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[i] Rainfall varies in space and time in a highly irregular manner and is described naturally 
in terms of a stochastic process. A characteristic feature of rainfall statistics is that they 
depend strongly on the space-time scales over which rain data are averaged. A spectral 
model of precipitation has been developed based on a stochastic differential equation of 
fractional order for the point rain rate, which allows a concise description of the second 
moment statistics of rain at any prescribed space-time averaging scale. The model is thus 
capable of providing a unified description of the statistics of both radar and rain gauge data. 

The underlying dynamical equation can be expressed in terms of space-time derivatives of 
fractional orders that are adjusted together with other model parameters to fit the data. The 
form of the resulting spectrum gives the model adequate flexibility to capture the subtle 
interplay between the spatial and temporal scales of variability of rain but strongly 
constrains the predicted statistical behavior as a function of the averaging length and time 
scales. We test the model with radar and gauge data collected contemporaneously at the 
NASA TRMM ground validation sites located near Melbourne, Florida and on the 
Kwajalein Atoll, Marshall Islands in the tropical Pacific. We estimate the parameters by 
tuning them to fit the second moment statistics of radar data at the smaller spatiotemporal 
scales. The model predictions are then found to fit the second moment statistics of the gauge 
data reasonably well at these scales without any further adjustment. 

Citation: Kundu, P. K., and J. E. Travis (2013), A stochastic fractional dynamics model of space-time variability of rain, 
J. Geophys . Res. Atmos., 118, 10,277-10,295, doi:10.1002/jgrd.50723. 


1. Introduction 

[ 2 ] Because of its irregular nature, rain is both difficult to 
measure accurately and predict from a physical model. 
Models of rain statistics provide a simple and conceptually 
economical way to capture the space-time variability of 
precipitation in terms of a small number of adjustable para- 
meters. They can be relatively easily validated from a large 
space-time data set, and once the parameters are tuned to 
data, the model provides a rather efficient method of des- 
cribing various statistical properties of precipitation over 
areas of similar rain climatologies. 

[ 3 ] In practice, rainfall is generally measured as a 
(nearly) instantaneous area-averaged quantity in radar 
measurements or as a time-averaged quantity at a point in 
rain gauge measurements. In theoretical models, they are 
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conveniently represented as suitable space and/or time ave- 
rages of a continuous stochastic field. This continuum appro- 
ximation is valid at the space-time resolution of the usual 
measurement methods under normal rainy conditions in 
which the inherent discreteness of rain at the scale of individ- 
ual drops is smoothed out. In this paper we develop a phe- 
nomenological model of space-time statistics of rain in 
terms of a random field R(x, t ) denoting the instantaneous 
point rain rate. It should be emphasized that R(x, t) is not 
directly observable, but when suitably area- or time-averaged, 
corresponds to measured quantities. 

[ 4 ] Radar scans yield area averages at an instant of time 
with horizontal spatial resolution of order 1 km. Rain gauge 
and disdrometer observations, on the other hand, lead to 
time-averaged point rain rate estimates with temporal 
resolution of order 1 min. These data can then be further 
“coarse-grained” by aggregating them to any desired larger 
space-time scale. Although rainfall varies in an apparently 
irregular manner, the underlying physical processes take 
place over an extended space-time region that causes the rain 
rate field to be correlated in space and time. Moreover, the 
statistics of rainfall depend on the space-time averaging 
scales in a nontrivial manner. In fact, it is well-known that 
[see, e.g., Bell , 1987] there is a subtle interplay between 
the space-time scales associated with the decay of the corre- 
lation function and the averaging scale that is reminiscent of 
fluid turbulence. Like the velocity field of a turbulent fluid, 
the rain rate field R(x, t) has the property that the larger the 
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averaging area, the longer the field remains temporally corre- 
lated. Similarly, the longer the period of time averaging, the 
greater the distance over which the spatial correlation per- 
sists. This property of the space-time correlation of rain is 
most easily captured in terms of the Fourier spectrum of the 
field. A spectral model of precipitation statistics was devel- 
oped in a number of earlier papers [Bell and Kundu , 1 996] 
(hereinafter BK96) and [Kundu and Bell , 2003] (hereinafter 
KB 03) that incorporates these features in a qualitative man- 
ner. The spectrum is generated from a Langevin-type stochas- 
tic differential equation for the spatial Fourier amplitudes 
of R(x , t) that is suggested by analogy with Brownian mo- 
tion. The model spectrum in turn directly determines the 
complete second moment statistics of the rain field averaged 
to any desired space-time scale and is thus in principle capa- 
ble of fitting both radar and rain gauge data. Thus, if the 
parameters of the model are tuned to fit the statistics of 
area-averaged rain rate using radar rainfall data, the model 
is then expected to describe, without any further adjustment , 
the statistics of time-averaged point rain rate data from a 
gauge network within the same space-time domain. We 
refer the reader to BK96 and KB03 for a more complete 
account of various aspects of the model and the relevant 
literature on other modeling approaches. 

[ 5 ] With the availability of large precipitation data sets in 
recent years, it has become feasible to validate the model 
quantitatively over a large range of space-time scales. Large 
multiyear data sets have now been produced from ongoing 
radar and rain gauge measurements collected as part of the 
ground validation (GV) program pursued by NASA during 
the Tropical Rainfall Measuring Mission (TRMM) [Wolff 
et al ., 2005]. In particular, a large amount of space-time 
colocated data is available from radar and gauge observations 
at several TRMM GV sites. We sought to test the spectral 
model described in BK96 and KB03 with the TRMM GV 
data. During this effort, it became clear that the model in 
its original form broadly captures the general features of 
the space-time statistics of radar-derived precipitation data 
but does not accurately fit the details. Moreover, if the 
model parameters are estimated by fitting the radar data, 
the predicted statistics of the accompanying rain gauge 
data depart substantially from the observed statistics. 
Alternatively, the parameters estimated independently from 
the radar and gauge observations belong to qualitatively 
distinct model regimes. The inevitable conclusion was that 
the originally proposed model spectrum needed to be gener- 
alized for it to describe both radar and gauge observations. 

[6] Our present work stems from an attempt to find such a 
generalization. With the integration of radar and gauge obser- 
vations, a larger range of space-time scales becomes experi- 
mentally accessible. In order to achieve greater flexibility in 
fitting all the available data, we extend the model framework 
by generalizing the underlying stochastic dynamical equation 
from an ordinary differential equation in time to a differential 
equation of a suitable noninteger order. The new model is 
able to fit the second moment statistics of both radar and 
gauge data more closely than the original model over the ac- 
cessible range of space-time scales. It should be emphasized 
that the model describes only the second moment statistics of 
R(x, t), not the full probability distribution which is also 
known to depend on the space-time averaging scale [Kedem 
and Chiu , 1987; Bell , 1987; Kundu and Siddani , 2007]. 


[ 7 ] The stochastic equation introduced in this paper to 
describe the precipitation process involves a mathematical 
framework generally referred to as fractional calculus. 
Broadly speaking, fractional calculus constitutes an exten- 
sion of the notion of derivatives and integrals of ordinary 
calculus to derivatives and integrals of fractional order [Miller 
and Ross , 1993; Oldham and Spanier , 2006; Samko et al ., 
1993]. West et al. [2003] has given a thought-provoking 
account of how such fractional operators can arise in the 
description of a wide variety of macroscopic physical pro- 
cesses. While these fractional differential operators can be 
mathematically formulated in several different ways, their 
representation, as certain integral operators with a power 
law kernel known as Riemann-Liouville operators, lends it- 
self to the clearest physical interpretation. The new model 
generalizes the “old” spectral model of BK96 and KB03 by 
replacing the ordinary time derivative in the underlying sto- 
chastic dynamical equation by a fractional derivative opera- 
tor. Because of the postulated power law dependence of the 
relaxation time of the Fourier modes in BK96 and KB03, the 
fluctuations of the rain field were already implicitly nonlocal 
in space. Now the nonlocality in time evolution implied by 
the power law kernel of the fractional time derivative operator 
also reflects the presence of a memory [Beran, 1994]. 

[8] One noticeable shortcoming of the original BK96 
model with the time evolution governed by a first-order time 
derivative was that the model did not fit the lagged autocorre- 
lation function of the area-averaged rain rate very well. The 
falloff rate of the lagged autocorrelation with lag z predicted 
by the model was found to differ markedly from what was 
actually observed, especially for small r, and it was suggested 
in BK96 that this indicated the need for introducing higher 
order autoregressive processes. The introduction of a frac- 
tional-order time derivative allows us to control the shape at 
small r effectively in a parsimonious manner and thereby ob- 
tain much better fit to the observed lagged autocorrelation in 
this regime. 

[ 9 ] We have tested the validity of our model using two radar 
data sets belonging to TRMM standard product 2A-53 that 
were generated as part of the TRMM GV program: (i) a 
spatially gridded set of images generated from scans by a 
National Weather Service radar (hereafter called MELB) 
located near Melbourne, Florida and also (ii) a similar data 
set from the radar (hereafter called KWAJ) located on the 
Kwajalein Atoll, Republic of Marshall Islands in the Pacific 
Ocean. The MELB radar has the advantage that the portion 
of its field of view (FOV) that is over land contains a dense 
network of rain gauges. However, its coastal location creates 
a somewhat complicated precipitation climatology. On the 
other hand, the KWAJ radar FOV has the advantage of being 
in a predominantly oceanic environment. However, the few 
gauges that are available in the area are rather sparsely distri- 
buted. The gauge data used in this paper are part of TRMM 
standard product 2A-56. 

[ 10 ] The remainder of the paper is organized as follows. In 
section 2, we give an account of the basic mathematical 
framework of the new stochastic model. In section 3, we first 
describe the radar data analysis. We then discuss the process 
of estimating the model parameters by fitting the model 
predictions to the observed second moment statistics of the 
MELB and KWAJ radar data. Finally, we test the model with 
the rain gauge data by examining how well the model tuned 
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to the radar data reproduces the second moment statistics 
of rain data derived from a cluster of gauges located within 
the radar FOV. Section 4 is devoted to a discussion of the re- 
sults along with the various caveats. The paper is concluded 
in section 5 with a summary of the findings and some direc- 
tions for future work. In Appendix A we give a brief account 
of fractional calculus. Appendix B presents some details of 
the mathematical derivations of the necessary formulas and 
will be frequently referred to in the main text. 

2. The Model Framework 

[ 1 1 ] The F ourier spectrum of the precipitation field provides 
a convenient way to characterize the various aspects of its 
space-time variability in a succinct manner. In this section 
we first constmct a stochastic dynamical model for the local 
rain rate field that naturally leads to such a spectrum. We then 
relate the space-time covariance statistics of the area- and time- 
averaged rain rate to the spectrum through the Fourier trans- 
form representation. The resulting formulas are derived in 
Appendix B. In the last two subsections, we examine the be- 
havior of the model in the limit of vanishingly small space- 
time scales. 

2.1. The Basic Equations 

[ 12 ] In this subsection we describe the basic theoretical 
framework of the spectral model. A central quantity of 
interest is the space-time covariance of the point rain rate 
field R(x, t) at points x, x f in a two-dimensional Euclidean 
plane (neglecting the Earth’s curvature) and at times t , t' 

c(x, t\ x', f) = (R'(x, t)R'(x', ?)), (1) 

where R\x, t)=R(x, t) — (R) is the deviation of the rain rate 
from the mean and the angle brackets (...) denote ensemble 
average over similar rain climatologies. In our model it is 
determined from the Fourier spectrum of the rain field. 

[ 13 ] As in BK96 and KB03, we assume the rain statistics to 
be spatially homogeneous, isotropic, and temporally station- 
ary, (for brevity collectively referred to as being space-time 
stationary). The homogeneity and stationarity assumptions im- 
ply that c(x, t\ x\t') depends only on the difference between 
the space and time arguments, i.e., the spatial separation vector 
p = x - x' and the lag z = t — f. Isotropy further restricts the 
dependence to the form 

c(x, x', f) = c(p , r), (2) 

where p = |p|. We should note that in the present paper the 
stationarity property refers only to the second moment statis- 
tics rather than the full underlying probability distribution. 
The spatial Fourier amplitudes 

a( k,t) = (27r)~ 1 fd 2 x e~ lk ' x R'(x, t) (3) 

are, in general, complex but constrained to satisfy the condi- 
tion a*(k,t)=a(—k,t) (where asterisk denotes complex 
conjugation), which follows from the fact that R'(x, t) is real. 
We assume that the a(k,t) evolve in time according to the 
generalized Langevin equation 

-00 D 1 «(k, t) = -r/a(k, t) +/(k, t ) . (4) 


Here _ x £>f denotes the Liouville-Weyl fractional derivative 
operator of order ft with respect to the argument t, which can 
be regarded as a shorthand for the operator defined in 
equation (A2) with the lower limit of integration tending to 
—oo. See Appendix A for a brief account of some basic results 
from the calculus of fractional derivatives. From the definition, 
it is clear that equation (4) is actually an integro-differential 
equation and therefore represents nonlocal time evolution. 
Also, in equation (4), / (k, t ) represents a white-noise forcing 
term with zero mean and ^-function covariance 

(/ (k, t)f*( k', t')) = (2*) 3/2 F 0 tf(k - k')d(r), (5) 

and 

T k = z 0 (l + ^L 2 0 y a/2 ( 6 ) 

is the relaxation time for the Fourier mode k depending only 
on the wave number k = |k| by virtue of spatial isotropy. [Note 
the incorrect normalization in KB03 equation (5).] In the 
frequency domain, equation (5) is equivalent to 

(/ (k, «)/*( k', co')) = (2^) 3/2 F 0 <5(k - k')S(co - co '), (7) 

/ (k, co) being the temporal Fourier transform of / (k, t). (We 
will occasionally denote a function and its Fourier transform 
by the same symbol when there is no possibility of confu- 
sion). Here F 0 is a strength parameter, and r 0 and L 0 are char- 
acteristic time and length scale parameters, respectively. The 
characteristic length scale L 0 effectively separates the rain field 
fluctuations into two regimes: a short wavelength (large k) 
scaling regime in which z k tends to zero according to a 
power law k~ a and a long wavelength (small k) regime in 
which z k approaches r 0 . Physically, r 0 represents the duration 
of an average rain event. Equations (4)-(6) are the basic 
equations of the model. The three quantities F 0 , r 0 , and L 0 
together with the two dimensionless exponents a and ft define 
the full set of model parameters. Selection of the lower limit of 
time integration as -oo is dictated by the condition that we want 
the model to describe stationary temporal statistics, where 
there is no preferred choice of the initial time. Choosing any 
finite value would correspond to choosing a particular time 
at which the initial condition for the stochastic equation has 
to be set, thus leading to loss of stationarity. The rain rate field 
R'(x , t) itself defined by the inverse spatial Fourier transform of 
(3) satisfies a stochastic field equation involving fractional 
spatial and temporal derivative operators. The spatial deriva- 
tive operator that results from a spatial Fourier transform of 
equation (4) with z k given by equation (6) can be formally 
represented as the familiar Helmholtz operator (— V 2 +Zq 2 ) 
raised to the power aft 12. 

[ 14 ] The “old” spectral model of BK96 and KB03 is recov- 
ered in the special case /? = 1 , in which the derivative operator 

t reduces to the ordinary time derivative d/d t. Equation 
(3) then simply becomes an ordinary first-order stochastic 
differential equation with exponential relaxation. For a gen- 
eral noninteger order, the power law kernel of the fractional 
derivative operator is indicative of an underlying random 
process that is non-Markovian. The process is now character- 
ized by a nonexponential relaxation. The response by the 
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precipitation process to a unit impulse is determined by the 
Green’s function G(k, t — t ') , which is the solution to the 
inhomogeneous equation 

(_oX + G(k, t-t')= d(t - ?). (8) 

[ 15 ] The Fourier amplitudes a(k, t) have zero mean and 
lagged covariance of the form 

(a( k, t ) a*(k', t')) = 2i zc(k, r)<5(k — k'), (9) 

where, as a consequence of spatial isotropy, we can write 
c(k,z) =c(k, r) for the spatial Fourier transform of c(p,z). 
In the frequency domain, equation (9) corresponds to 

(a( k, co) a*( k', of)) = (2i z) 3 ^ 2 S(k, a>)J(k — k')<5(co — cot), (10) 


[ 17 ] The space-time covariance c(p, r) is given by the spa- 
tial Fourier transform of c(k,z). Spatial isotropy allows one 
to carry out the angular integration in the k-plane. This yields 

c(p, r) = Jo d£ kJ 0 (kp)c(k , z ), (14) 

where Jo(x) is the usual Bessel function of order zero. The 
analytical form (13) implies that the space-time covariance 
function c(p,z) given by (14) is in general not factorizable 
into spatial and temporal dependence. Upon setting r = 0, 
we obtain the spatial covariance function c(p , 0) in the form 
of an integral that is identical with the one encountered in 
the = 1 case. It can be similarly evaluated with the result 

c(p,0) = y 0 Cv(p/L 0 ), (15) 

provided we now identify the index v through the formula 


where we have introduced the power spectrum of the rain 
field fluctuations S(k,co) as the temporal Fourier transform 
of c(k,z). From the correspondence _ 00 Z)fo (— ioof under 
the action of the Fourier transform (see Appendix A), it 
follows that 


«( 2/J-l) = 2(l+v), (16) 

where, as in BK96 and KB03, we introduce 

C v (z) = (z/2) v K v (z), (17) 


S(k , co) 



2 ’ 


which can be simplified to the form 

S(k, co) = F 0 \\co^ + 2 cos(/?;r/2)|ew|^r^ + . (11) 

The principal branch of the multivalued function in the denom- 
inator is assumed to be in the range —n< Arg co<n. 

[ 1 6] The spectrum S(k, co) given by equation (11) yields the 
full set of second moment statistics of the rain rate field. The 
point covariance function c(p,z) is the space-time Fourier 
transform of the spectrum, i.e., 


c{p,z) = (2n) 3 ^ 2 J_ 00 dft;Jd 2 k e l ^ p mz>> S(k,co). (12) 


K v (z) being the modified Bessel function of order v. For a 
purely spatial random process, a class of covariance func- 
tions of the form (15) was originally introduced by Matern 
[1960,1986]. The factor y 0 is now a slightly more compli- 
cated quantity with the dimension of [rain rate] 2 that is 
expressible in terms of the basic model parameters: 


L 2 0 r(i+v) ’ 


where T(z) denotes the Euler gamma function. The temporal 
covariance of the rain rate field can be obtained by simply 
setting p = 0 in equation (14): 


fo° 

c(0, z) = J 0 dk kc(k , z). 


It can be evaluated in two steps. First, consider the temporal 
Fourier transform of S(k, co): 


c{k , z) = (2n) 



d 00 e ia>z S(k 1 co). 


From the explicit form (11) of the spectrum, it follows, 
by a simple scaling argument, that c{k,z) has the func- 
tional form 


c(k,t) =g(P)F 0 tf x h(\x\/x k ), (13) 

where h(x) is a certain transcendental function defined in 
Appendix B and g(fi) is a normalization factor adjusted 
so that h( 0)=1. When /?=1, h(x) is simply exp(— x). 
Unfortunately in the /? 7 ^ 1 case, it does not appear possible 
to express h(x) in terms of familiar analytical functions. 
The forms of the function h(x) for some typical values 
of P obtained by numerical integration are shown in 
Appendix B (see Figure Bl). 


[is] The point variance c(0, 0) = o\ is evaluated by letting 
p— >0, z — > 0 in equation (14) and making use of equation 
(13). In terms of the dimensionless variable y = 1 + £ 2 Zq, 
it takes the form 


^0 


= (i/2)y 0 r(i + 


f 00 

v)Ji d yy 


(l+v) 


(19) 


The integral converges when v > 0, i.e., when a(2fi — 1) > 2 
yielding o\ = y 0 r(v)/2, but diverges when v<0 causing 
ctq to be infinite. As was already discussed in BK96 and 
KB03, when v < 0, the spatial covariance c{pfi) has a power 
law singularity atp = 0: c(p,0) ~ (p/L 0 )~ which weakens 
to logarithmic when v = 0 [North and Nakamoto , 1989]. 
Later in section 3, as we examine the statistical properties 
of precipitation data sets, it will become clear that the model 
fit to data strongly favors the v < 0 case. 

[ 19 ] While the point rain statistics themselves cannot be 
directly measured, they determine the covariance statistics 
of the space- or time-averaged rain rate data that can be mea- 
sured. The next two subsections summarize their properties. 
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2.2. Covariance Statistics of Area-Averaged Rain Rate 

[ 20 ] Now we consider the space-time covariance of spa- 
tially averaged radar rainfall data. A gridded radar precipita- 
tion data set consists of a sequence of images, each image 
being an array of L x L square grid boxes in which the aver- 
age (near-instantaneous) rain rate is specified. Let 

r A {t)= (1/I 2 )jd 2 xi?(x,/) (20) 

A 

denote the average rain rate in an L x L square A at time t. The 
space-time covariance between the rain rates in two squares A 
and A' of equal area L 2 , whose centers are separated by a dis- 
tance vector s at two different times t and t+ z, is defined as 

r A4'( S > T ) = { r 'A(ty A '(t + T)), (21) 

where, as before, the prime on the rain rate variables indicate 
deviation from the mean. It can be written as a double area- 
integral over the space-time covariance function of the point 
rain rate: 

r A4 '( s > T ) = (!/ i4 )l d 2 x 1 d 2 x' c( s + X' - x, r). (22) 

^ A' 

Evaluating equation (22) for s = 0, r = 0, one obtains the var- 
iance of area-averaged rain rate cPym (r' A ). The spatial corre- 
lation between two boxes A and A’ separated by s is then 
given by O^'(s) = r^(s, 0)/<r 2 . 

[ 21 ] The time dependence of the lagged autocorrelation 
function ® aa (t) = Taa (0, r) / a A defines a characteristic time- 
scale, the integral correlation time for area-averaged rain rate 

t a = J 0 dr 0^4 (r). (23) 

For an exponentially decaying correlation, z A is simply the 
(1 /^-folding time. However, our model-predicted correlation 
functions are markedly nonexponential. 

[ 22 ] Explicit formulas for the various second moment sta- 
tistics including a A and t a are given in Appendix B. 

2.3. Covariance Statistics of Time- Aver aged Rain Rate 

[ 23 ] Next we examine how the model represents the sta- 
tistics of rain gauge data. While the actual quantity mea- 
sured depends on the specific type of instrument, the data 
is usually converted into a form that can be idealized as 
the time-averaged rain rate at a point, i.e., 

r r (x) = (l/r)J o r d/ J R(x,?). (24) 

In general, one could consider the space-time covariance Y TT ' 
(p, r) between rain rates at two points separated by a distance 
p, time-averaged over two intervals T and T whose mid- 
points are separated by a lag r (see Appendix B). For conve- 
nience we restrict ourselves to the zero-lag case, namely the 
spatial covariance of the rain rate averaged over a time inter- 
val [0,7] at two points with separation p 

Trr(p,0) = (rV(x)rV(x')) 

= (i/r 2 )J 0 dJ 0 d ? ' c(p,t-f), 


which can be reduced to a single integral 

rirfoO) = (2/nf o r dT(l - t/T) c(p,t). (25) 

The limit p — > 0 yields the variance of time-averaged rain rate 

= (2/n C dr(l - t/T) c(0, r). (26) 

[ 24 ] From equations (25) and (26) we can compute the spa- 
tial correlation function of gauge pairs as a function of sepa- 
ration: x Y T t{p) = Tjt(p, 0)/g 2 t . For explicit evaluation of 
these quantities, it is again convenient to go over to the 
Fourier representation. The resulting formulas are given in 
Appendix B. 

2.4. The Power Law Scaling Regime in the Case v < 0 

[ 25 ] As already noted, the space-time covariance c(p,z) 
predicted by the model becomes singular in the limit of small 
p and r when v < 0 but approaches a well-defined finite value 
when v > 0. In the case v < 0, which is of greater interest to 
us, this singular behavior at small space-time separations p, 
r indicates the presence of a “scaling regime” in which the 
model exhibits approximate scale-invariance under a space- 
time scaling p — ► Xp, z — > X a z. This is elucidated by examining 
the limiting form of the spectrum S(k , co) in the limit of large k 
and co, which we denote by co). We find that S^°\k, co) is 
invariant (up to an overall multiplicative factor) under a scale 
transformation in the Fourier space k^X x k, co ~^X~ a co: 

S (oo) [X~ x k, X~ a co) = X 2aP S^ (k, co). (27) 

Its Fourier transform c (oo \p,r), which is the asymptotic form 
of the exact space-time covariance c(p,r), satisfies 

c^\Xp,X a z) = r 2|v| c^(p,r). (28) 

As the scale factor 2 — > 0, one attains larger and larger values 
of k and co corresponding to smaller and smaller space-time 
scales. Choosing the scale factor to be X = 1/p*, we immedi- 
ately conclude from equation (28) that c (oo) (p,r) must have 
the functional form 

c (oo) (p,r) = y 0 p^~ 2lvl (p(n/p a ^aj), (29) 

where, for convenience, we have introduced the dimension- 
less variables p*=p/Z 0 , = r/r 0 . The type of combined 

space-time scaling property expressed by equation (28) is 
frequently referred to as dynamic scaling , a being the corre- 
sponding scaling exponent. The scaling function (p(£\ a , jJ) 
also depends explicitly on both a and /?. The limiting 
behavior of c^(p, r) as p*, r* — > 0 is nonuniform, i.e., 
dependent on the direction from which the space-time 
origin is approached. The double limit is characterized by 
the functional dependence of cp{^; a , ft) on the scaling vari- 
able £ = r^/p^. In particular, the behavior of c^(p,z = 0) 
asp* — > 0 is determined by the asymptotic behavior of cp(^; a, ft) 
near £ = 0. Similarly, the behavior of c ((X) \p = 0, r) as r* — > 0 is 
determined by the asymptotic behavior of (p(g; a, fX) as 
£ — > 00 . From the exact result for the spatial covariance 
c{p,z = X)) = yQC v (p*) [see equation (15)], making use of 
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Figure 1. Map showing the locations of the MELB radar and the tipping bucket rain gauges in Florida 
operated as part of the TRMM ground validation program. Also shown are the radial distance contours 
of the radar FOV and the position of the 128 km area from which statistics is collected. 


the asymptotic behavior of C v (z) as z— >0 when v<0, we 
get the limit 

(p(0; a,fi) =2-( 1+2 l v Dr(|v|). (30) 

[ 26 ] Unlike the /? = 1 case studied in BK96 and KB03, 
c(p = 0, r) likely cannot in general be expressed in closed form. 
However, scaling arguments like the one above can also be in- 
voked to obtain the asymptotic r-dependence of c^\p = 0, r) as 

r — > 0. Setting p = 0 and choosing the scale factor X = x~^ a in 
equation (28), we get the asymptotic form 

c M(p = o,6~yo^“ 2|v|/a (3i) 

as r*— >0, where K is a dimensionless constant. This is 
compatible with the general form of equation (29) if and only 
if (p(£\ a, /?) has the falloff 

p(£;a,jS) ^ J &r 2|v|/a . (32) 

f->°o 

The asymptotic behavior of the area and time-averaged 
statistics can be deduced by replacing the exact space-time 


correlation c(p,z) by its scaling approximation c (oo) (p, r) in the 
area/time integrals. The details are relegated to Appendix B. 

[ 27 ] The scaling properties of the space-time covariance in 
the p=\ case were explored in greater detail by Kundu and 
Bell [2006]. 

2.5. Effect of a Short Distance Cut-Off 

[ 28 ] As will be found later, a problem arises when we 
examine the spectral model predictions in light of the data. 
We find that the growth property of the space-time covari- 
ance at small space-time scales in the v < 0 case predicted 
by the spectral model is broadly consistent with data, but 
only up to a certain point. The spatial variance o 2 A estimated 
from radar data does appear to follow the predicted power 
law behavior [see equation (B22)] within the accessible 
range of spatial scales. But rain gauge data, which allows 
one to access much smaller space-time scales, exhibit a more 
complicated behavior. The scale dependence of the temporal 
variance o\ follows the model prediction up to a certain 
minimum averaging time T of about 30 min, but increasingly 
deviates from it at shorter time scales. As T 7 — ► 0, instead of 
growing in accordance with the predicted power law growth, 
a? appears to gradually approach a finite asymptotic value o\, 
the point variance, indicating a steeper decrease of the 
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Table 1. Model Parameters and Related Quantities for the Spectral Model 


Radar 

Season 

(S) (mmh -1 ) 

a 

p 

V 

Jo (mm 2 h 2 ) 

T 0 (km) 

r 0 (min) 

cr 0 2 (mm 2 h 2 ) 

A (km) 

KWAJ 

MAM 2001 

0.098 

0.99 

1.18 

-0.327 

0.019 

281 

775 

2.5 

0.48 

KWAJ 

JJA 2001 

0.232 

0.93 

1.28 

-0.279 

0.060 

438 

770 

7 

0.36 

KWAJ 

SON 2001 

0.357 

0.99 

1.24 

-0.265 

0.213 

136 

381 

13 

0.27 

KWAJ 

DJF 2002 

0.114 

1.40 

1.00 

-0.298 

0.067 

72.1 

524 

3.5 

0.32 

MELB 

DJF 2001 

0.028 

1.17 

1.18 

-0.202 

0.030 

69.0 

196 

1.2 

0.07 

MELB 

MAM 2001 

0.104 

1.17 

1.26 

-0.113 

0.348 

73.2 

197 

6 

0.09 

MELB 

JJA 2001 

0.214 

1.14 

1.26 

-0.130 

1.078 

33.9 

98.8 

13 

0.19 

MELB 

SON 2001 

0.183 

1.12 

1.20 

-0.218 

0.337 

51.5 

209 

10 

0.18 


characteristic time scale i k than the model-predicted k~ a 
falloff at large k. On the other hand, the v > 0 case of the 
model, which would have accounted for this effect, does 
not generally fit the radar statistics. A simple way to incor- 
porate the desired modification at short length and time 
scales is to assume that instead of equation (6), T k is given 
by the relation 


n 


z 0 (l+k 2 L 2 0 ) a/2 ;k<l/A 
0 ;k>l/A 


(33) 


to obtain an estimate of A, we return to equation (19) for the 
point variance. Explicit evaluation of the integral with 
equation (33) leads to the simple analytic result: 

=^ol r WI'[(l + L 2 /A 2 ) W - l], (34) 

The spatial cut-off A is easily obtained from equation (34) if a 
reasonable value of o\ can be estimated by extrapolation from 
the gauge data. 


where A is a short distance (“ultraviolet”) cut-off. Physically 
this means that the spatial Fourier modes of the precipitation 
field a(k, t) of wavelength shorter than In A are damped out 
instantly. Accordingly, the spectrum S(k , co) and its tempo- 
ral Fourier transform c(k,r) also vanish for wave numbers 
k> HA. In a more realistic model that describes the small- 
scale behavior of rain, the sharp cut-off introduced in 
equation (33) may have to be appropriately modified based 
on empirical evidence. The possibility of a steep decrease 
of the spectrum beyond the scaling regime is reminiscent 
of a similar phenomenon in fluid turbulence, namely break- 
down of the famous Kolmogorov scaling of the energy 
spectrum beyond the inertial range [Frisch, 1995]. In order 


3. Comparison Between the Model 
and Observations 

3.1. Radar Data Analysis 

[ 29 ] The parameters of the model can be conveniently esti- 
mated from a gridded radar precipitation data set. We have fit 
the model to the data sets (TRMM standard product 2A-53) 
constructed from radar images collected from the KWAJ and 
MELB radars located respectively at (8.718°N, 167.733°E) 
and (28.1 13°N, 80.654°W) (Figure 1). Next we summarize 
some relevant features of these data sets. 

[ 30 ] For both the KWAJ and MELB radars, the FOV con- 
sists of a circular area of diameter 300 km. The data was 
gridded at a 2 km x 2 km spatial resolution. In order to reduce 







Sox size L (km) 


Figure 2. Plot of the variance of area average rain rate <7 a - (^(L) as a function of the spatial 
averaging scale L estimated from (left) the Melbourne radar data for the four seasons Winter 
(December 2000-February 2001), Spring (March-May 2001), Summer (June-August 2001), and 
Autumn (September-November 2001) and from (right) the KWAJ radar for the four seasons 
Spring (March-May 2001), Summer (June- August 2001) and Autumn (September-November 2001), 
and Winter (December 2001-February 2002) superimposed on the curves predicted from the model 
formula (B7) with the parameters listed in Table 1. 
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Figure 3. Comparison of the spatial correlations of 2 km radar pixels as a function of the separation s as 
estimated from the Melbourne radar data (the + symbols) and the spectral model (the solid curve) with the 
parameters listed in Table 1 for four seasons spanning the period December 2000 to November 2001 as in 
Figure 2. 


the uncertainties in rain retrieval due to radar attenuation with 
distance, statistics were computed from a 128 km x 128 km 
area centered at the radar location by aggregating the data 
at spatial scales L = 2, 4, 8, ... , 128km. Only those grid 
boxes within the selected area that had at least 95% valid 
pixels were included in the sample. This helped eliminate 
boxes at the smaller scales 4, 8, and 16 km located near the 
center, which occasionally suffered from data dropout be- 
cause of ground clutter. 

[ 31 ] The temporal sampling pattern of the radar has an im- 
portant effect on the estimation of rain statistics. The MELB 
radar scans were carried out at different speeds during the 
quiet and the active periods of rainfall. Rainfall images were 
available from it every 5 min during the active periods and ev- 
ery 10 min during the quiet periods. Images from the KWAJ 
radar, which was operated at a single speed, are spaced in 
a somewhat more irregular manner: the temporal spacing 
between the consecutive images predominantly alternated 
between 5 and 7 min long gaps. 

[ 32 ] The various space-time statistics of the TRMM 2A-53 
data needed for the analysis (the variance and covariance 
functions) were computed using routines provided in the R 
statistical software package. For the spatial analysis we com- 
puted the variances o\ for various box sizes L between 2 and 
128km and the spatial autocorrelation function <D^'(s) for 


pairs of Z = 2km radar pixels. For simplicity the measured 
spatial correlation function is regarded as a function of the 
separation = |s| in accordance with our spatial isotropy as- 
sumption. Since all the available pairs within the 128 km area 
are included in the computation of O^'(s), our data analysis 
automatically averages over all pair locations within the box 
and allows us to assume spatial homogeneity. The direction 
dependence was implicitly averaged over two mutually per- 
pendicular directions in this step by pooling together all pairs 
with the same separation s. 

[ 33 ] One effect of the uneven temporal sampling of different 
portions of a data set like that of the MELB radar is that it in- 
troduces a systematic bias into the usual estimates of many sta- 
tistics. In particular, the “simple” estimates of the elementary 
statistics at a spatial scale L , such as the unconditional mean 
( R ) and the probability of nonzero rain p{L) = Pr[r A > 0] 
obtained by summing the corresponding data and dividing 
by the actual number of samples N(L ), are biased high. One 
can obtain “improved” estimates of these statistics that reduce 
this bias as follows. We assume that the temporal spacing 
between two successive images At 2 during a quiet period 
(labeled 2), which is mostly dry, is larger than the spacing 
At! during an active period (labeled 1) when rain occurs more 
frequently. Suppose At 2 = wAt\; in the case of the MELB 
radar, w = 2. For each stretch of the time series of radar images 


10,284 


KUNDU AND TRAVIS: A SPECTRAL MODEL OF RAIN STATISTICS 


3 , 

< 


e 


c 

g 

M. 

o> 


o 


O 

aj 


CO 



Separation s (km) 


Figure 4. Same as in Figure 3 but for the KWAJ radar data (the + symbols) and the spectral model (the 
solid curve) for the four seasons spanning the period March 2001 to February 2002. 


during a quiet period, we augment it by padding each gap 
successively with a segment consisting of (w — 1) additional 
copies of the initial point spaced so that the new series has a 
constant time step At\. 

[ 34 ] We can get a rough estimate of the bias from uneven 
sampling in the following way. Let us write N=N\ +N 2 , 
where N\ and N 2 are the actual number of samples during 
the active and the quiet periods, respectively. For the nonzero 
rain probability p , instead of the naive estimate 


Aj >0 + A ^ 0 

N\ +N 2 


( 35 ) 


we now have the adjusted estimate 


(36) 

where A^ >0 andA^ 0 are the number of rainy samples during 
the active and quiet periods, respectively, and N'=N\ + wN 2 
is the total number of samples in the augmented series. If 
the number of rainy samples occurring during the quiet 
period is small, i.e., N t 2 >0 «N\ then we get the approximate 
result p' ~ (N/N r )p. The same correction factor (AW) approx- 
imately applies to the unconditional mean (R) and the vari- 
ous unconditional moments M q (L) = (rf) under the similar 
condition that the contribution from rain occurring during 
the quiet periods can be neglected. 


[ 35 ] The lagged autocorrelation O aa ( z) was computed for 
several averaging scales L between 2 and 128 km and for var- 
ious lags t from the sequence of radar images. As expected, it 
was found that observed autocorrelation in general exhibits 
fluctuation due to sampling uncertainty, which becomes large 
especially at those r values where one has relatively few sam- 
ples. For our purpose, however, only those lags for which 
the sampling error is small enough to be negligible were 
considered. In order to reduce the effect of the sampling er- 
ror, we had to restrict ourselves to lags for which the number 
of available samples exceeds a judiciously chosen value. 
The estimates at the smaller scales L = 2, 4, 8,. . ., 64 km were 
obtained by partitioning the 128 km area into nonoverlapping 
subareas of size L*L and averaging over the computed lagged 
autocorrelation for each subarea. Pooling together samples of 
area-pairs for a fixed lag r for the temporal autocovariance also 
presumes temporal stationarity of the statistics. 

3.2. Fitting the Model to Radar Rainfall Statistics 

[ 36 ] Estimation of the model parameters presented a diffi- 
cult problem even in the simpler f}= 1 case considered in 
BK96. Not surprisingly, the general /? case proves to be even 
more challenging. Since the parameters appear directly in the 
spectrum, one’s first thought might be to perform a spectral 
analysis of the pixel data and estimate the spectrum of the 
point rain rate from it. One could then fit the observed spec- 
trum to equation (11) and estimate the model parameters. 
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Figure 5. Comparison of the lagged autocorrelation at four different spatial scales (L = 2, 8, 32, and 
128 km) as a function of the separation s as estimated from the MELB radar data (the symbols) and the 
spectral model (the solid curve) with the parameters listed in Table 1 for the same period as in Figure 3. 


However, as was noted in BK96, this approach has the disad- 
vantage that the spectral estimates often suffer from distor- 
tions at spatial scales comparable to the size of the area 
considered (L= 128 km in our case) that are artifacts of the 
periodic boundary condition imposed by the analysis method- 
ology. This is an important issue for us since an intended appli- 
cation of the model is to perform a statistical intercomparison 
of rain rate data from ground radar and a rain gauge network 
located within the radar FOV. While there are smoothing tech- 
niques for the numerically evaluated spectrum that reduce 
these distortions somewhat [Jenkins and Watts , 1968], it is 
not clear to what extent the choice of a particular smoothing 
prescription would affect our results. Moreover, the radar data 
was not uniformly spaced in time, which would therefore ne- 
cessitate further interpolation. To circumvent these difficulties 
we choose to follow the approach adopted in BK96 and deter- 
mine the parameters by directly fitting the statistics of the area 
averages to the model predictions. Our fitting process is in 
effect suitably weighted so that the model faithfully represents 
the observed statistics at spatiotemporal scales of interest to us. 

[ 37 ] Our parameter estimation method takes advantage of 
the mathematical structure of the model and proceeds in 
two stages. First, the parameters y 0 , L 0 , and v are estimated 
from the spatial statistics. The parameters v and L 0 are 
obtained from a suitably weighted least squares fit to the spa- 
tial correlation function ® p i x (s) = (s) for 2 km pixels sep- 

arated by a distance s = |s|. For simplicity, the slight direction 


dependence of ® p i x (s) expected from the model at small sep- 
arations is neglected. From the gridded radar images, we first 
compute the (Pearson) sample correlation coefficients p m for 
all pairs of pixels with spatial separation vectors s m of length 
s m ( m= 1,2, ... , M) as the estimate of the theoretical spatial 
correlation ® p i x (s m ) for a set of M= 66 selected values of s m 
ranging from the minimum value of 2 km up to about 1 60 km 
(slightly less than the length of the diagonal of a 128 km box). 
The vectors s m are chosen so that all directions are more or 
less uniformly represented and there are a substantial number 
(in the range 200-23,000) of pairs N m . We then seek to obtain 
estimates of v and L 0 by minimizing the quantity 

J(yM) ^j m =l Wm \Pm — ®pix( s ra)] 5 (32) 

where w m are a set of weights inversely proportional to the 
variances 27 m 2 of the sample spatial correlation p m . If the sam- 
ples entering into the computation of p m are independent, then 
27 m 2 oc 1 !N m and consequently the weights are expected to be 
proportional to the number of samples, i.e., w m 00 N m . In real- 
ity, there is data dependency caused by the space-time corre- 
lation of the rain field, thus reducing the effective number of 
independent samples. The overall effect of this dependency 
for our estimation problem is somewhat difficult to assess 
and will not be taken into consideration here. Evaluation of 
the model spatial correlations ® p i x (s m ) is carried out from 
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Figure 6. Same as in Figure 5 but for the KWAJ radar data (the symbols) and the spectral model (the solid 
curve) for the same period as in Figure 4. 


equation (B5) using the software Mathematica, version 9, 
from Wolfram, Inc. The minimum of J(v, L 0 ) is sought by 
employing the numerical implementation of the Nelder- 
Mead Simplex Algorithm [Nelder and Mead , 1965; Press 
et al , 1992; Lagarias et al , 1998] in Mathematica. A starting 
value for the index v can be obtained by fitting the asymp- 
totic form equation (B22) to the variance a A for small L. 
The computation of the spatial parameters v and L 0 for each 
data set was quite efficient and took only a few minutes on 
a desktop computer equipped with two Quad-core Intel 
Xeon processors. The normalization constant y 0 is fixed by 
a least squares fit of the model prediction (using the exact for- 
mula (B7)) to the observed variances a A . The Nelder-Mead 
Algorithm, despite its heuristic nature, is our preferred choice 
for an optimization method primarily because of its simplic- 
ity and speed of execution. For a test run on a KWAJ radar 
data set, the Differential Evolution Method [Storn and 
Price , 1997], a global optimization method available within 
Mathematica, gave results for v and L () that are substantially 
identical with those from the Nelder-Mead Method but took 
much longer (about 28 times) to converge. 

[ 38 ] The next step is to estimate the temporal parameters 
r 0 , ft (and therefore a). They can be determined by fitting 
the lagged autocorrelation function Q>aa(j) for a particular 
box size L. We calculated the model values ® aa(?j) using 
equations (B6) and (B7) and fit them to the observed values 
Pj for a 1 6 km box and a set of time lags t = Zj (j = 1 , 2, . . . , n) 


up to a certain reasonably large r. Values of ft and r 0 are 
estimated by minimizing the quantity 

J{p, to) = [fij - (tj)] 2 ( 38 ) 

using the Nelder-Mead Simplex Algorithm in Mathematica. 
The weights Wj can now be taken to be 1 (the ordinary least 
squares case) since the numbers of samples are roughly equal 
for all lags in the range considered. The observed lagged 
autocorrelation was found to exhibit large systematic depar- 
ture from the model-predicted form for larger values of z. 
We therefore had to restrict the range of r over which the 
fit was carried out (about 400 min for KWAJ and 200 min 
for MELB). 

[ 39 ] Unfortunately, evaluation of the quantities ® 'aa(?J) in- 
volves computing (2+1) dimensional integrals in the Fourier 
space of k and co making the optimization problem computa- 
tionally rather onerous. The Nelder-Mead Algorithm took a 
long time to converge. The execution time depended on the 
number of lag values employed to define the sum of squares 
J (/?, ro); for each data set it typically took about 1 h per value 
of j on a workstation equipped with an Intel i7 870 processor. 
We should note that in estimating these parameters, the short 
distance cut-off parameter A is implicitly set to zero with the 
anticipation that it is sufficiently small compared to the 
smallest spatial scale accessible to radar measurements. Its 
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Figure 7. Plot of the variance of time-averaged point rain rate a T 2 as a function of the averaging scale T 
for precipitation data from the TB gauges located within the MELB radar FOV. The panels exhibit results 
for the same time period December 2000 to November 2001 divided into four seasons as in Figure 2. The 
open circles denote the observed variance and the solid curve represents the model results. 


value can only be meaningfully estimated from the gauge 
data at very short time scales. 

[ 40 ] We divided the 2A-53 annual radar precipitation data 
sets from the KWAJ and MELB radars into four 3 month 
long seasons — Winter (December-February, or DJF), Spring 
(March-May, or MAM), Summer (June-August, or JJA) 
and Autumn (September-November, or SON), and carried 
out the estimation of the model parameters for each. In doing 
so, we implicitly regard the statistics of each season as station- 
ary in the wide sense. Here we present results of parameter 
estimation for four seasons for each of these two radars. For 
the MELB radar, we selected the 2001 DJF, MAM, JJA, and 
SON seasons. For the KWAJ radar, DJF 2001 data had to be 
excluded from the analysis because of the highly erratic behav- 
ior of the temporal correlation, and so we considered DJF 2002 
instead as a representative example. 

[ 41 ] The model parameters are listed in Table 1 . The spatial 
index v lies in the range —0.27 to —0.33 for the KWAJ radar 
and —0. 1 1 to —0.22 for the MELB radar. The temporal index 
P is found to be substantially greater than 1 , in fact lying in 
the narrow range 1 .2-1.3, in all cases except one, namely 
the KWAJ DJF 2002 season. The characteristic length and 
time scales L 0 and r 0 are generally longer for the KWAJ data 
sets than the MELB data sets, which is presumably attribut- 
able to the purely oceanic environment of the former. The 


constant y 0 varies widely from one data set to another and ap- 
pears to be related to the mean rain rate. Figure 2 shows a 
comparison of the variances a A deduced from the data (sym- 
bols) and those computed from the model (solid curve) for 
various box sizes L. The fits to the spatial correlation function 
O aa (s) for the L = 2 km radar pixels as a function of the sep- 
aration s are shown in Figures 3 and 4 for the MELB and 
KWAJ radars, respectively. The quality of the fits for the 
temporal statistics is illustrated in Figure 5 (for the MELB 
radar) and Figure 6 (for the KWAJ radar) by the plots of 
® aa( r) at several other spatial scales. It is seen that the model 
fits the observation reasonably well for spatial scales L 
between 2 and 32 km, but the fit becomes worse at larger 
scales. By allowing the order of the time derivative ft to be 
an adjustable noninteger parameter, a closer fit at small z is 
achieved for the temporal statistics compared to the /? = 1 case 
of the model, as illustrated in Appendix B (Figure B2) for the 
KWAJ JJA 2001 season. 

3.3. Validating the Model with Rain Gauge Data 

[ 42 ] In this section we compare the model predictions 
for the statistics of time-averaged point rain rate outlined 
in section 2.4 with the statistics of rain gauge observa- 
tions. The temporal statistics were obtained from a rain 
gauge data set (TRMM standard product 2A-56) for the 
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gauge pair separation p {km) 

Figure 8. Plot of the spatial correlation ^tAp) of the daily averaged gauge data as a function of the gauge 
separation p for the same data set as in Figure 7. The panels exhibit results for the same time period 
December 2000 to November 2001 divided into four seasons as in Figure 2. The + symbols denote the ob- 
served correlation for each gauge pair. The open circles denote the observed correlation averaged over all 
gauge pairs in each distance bin interpolated by the solid curve. The error bars are estimated according to 
the method described in the text. The dashed curve represents the results computed from the model. 


same time period over which the spectral model parame- 
ters were obtained in section 3. The data consists of 
estimates of 1 min rain rate derived from rainfall accumula- 
tions in the tipping bucket (TB) gauges located within the 
same central 128 km area within the radar field-of-view 
(FOV) that was used to extract the radar statistics during a 
particular period of observation (Figure 1). For the KWAJ 
radar FOV, there were only seven gauge locations, each 
equipped with paired TB gauges and, moreover, not all of 
them were available for every season. For the MELB radar, 
there were about 50-60 gauges available within the central 
128 km box. 

[ 43 ] We computed an estimate of the observed temporal 
variance g 2 t for T in the range 1-7200 min by taking the 
average over a moving window of length T for each gauge 
along the time series and computing the variance for the full 
set of gauges for the entire season. The variances for the indi- 
vidual gauges were averaged together to estimate g 2 t . The 
spatial correlation of the gauge pairs ^tAp) was also com- 
puted for daily averaged rain rates (T= 1440 min). For the 
MELB data set, the average spatial correlations between 
gauge pairs were estimated by binning them into spatial 


separation bins. The average over each bin was estimated 
as follows. A Fisher z-transform z^tanh -1 ^, which maps 
the interval (—1, 1) onto the real line, was applied to the cor- 
relation estimates for each gauge pair. It is well-known that 
l Hawkins , 1989] for a fairly large class of probability distri- 
butions of the correlation estimates x ¥, the variable z is nearly 
normal and consequently, the standard estimates of the 
confidence intervals apply. The inverse transform x ¥ = tanh z 
is then applied to the statistics obtained from the normal 
theory to determine the average correlation and the corre- 
sponding 95% confidence intervals for each separation 
bin. In practice, the histograms of the z-transformed values 
turned out to be slightly non-normal. We therefore checked 
the accuracy of our results by constructing a bootstrap 
distribution from 10,000 replicas of a random sample drawn 
from the underlying empirical distribution and constructing 
bootstrap confidence intervals [Efron, 1981]. The resulting 
bias of the bootstrap estimate of the mean was found to be 
negligibly small and the estimated confidence intervals 
agreed well with the normal statistics estimates. For the 
KWAJ data, the gauge distribution was very sparse, and with 
the small number of rain gauges that were available in the area, 
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Figure 9. Plot of the (left) variance ot as a function of the averaging scale T and the (right) spatial cor- 
relation A* T7 (p) of the daily averaged gauge data as a function of the gauge separation p for the Summer 
(June- August) 2001 season KWAJ gauge data. 


only a few appeared to give reliable results. As a result, each 
pair was considered individually even though the resulting spa- 
tial correlation estimates were much noisier. Moreover, it was 
decided to leave out of the analysis certain gauges that showed 
anomalously small cross correlation with their duplicates and 
other nearby gauges. 

[ 44 ] It should be noted that the TB estimates at small T 
(less than about 7-10 min) are known to be somewhat 
unreliable [Habib et al ., 2001]. This is because the TB gauge 
data consist of a sequence of tip times of a bucket of a certain 
fixed capacity. The interpolation scheme based on a cubic 
spline algorithm [ Wang et al., 2008] used to generate the 1 
min time series of rain rates, leads to approximations when 
isolated single and double tips occur at low rain rates. The 
statistics of rain rates at small time scales T can thus be 
potentially distorted by the artifacts of this algorithm. 

[ 45 ] Using the parameters extracted from the fit to the 2A-53 
radar data, we compute the statistics of time-averaged point 
rain rate from the model using equations (B18) and (B19) 
and compare them with the corresponding observed quantities 
computed from the 2A-56 rain gauge data. This comparison 
provides a powerful validation of the model. Figure 7 exhibits 
our results for the observed and the model-predicted variance 
g 2 t as a function of T for the MELB DJF-SON 2001 data in 
the range 1-7200 min. A comparison of the corresponding 
spatial autocorrelation ¥ T7 {p ) of the daily averages as a 


function of gauge separation p for the MELB data is shown 
in Figure 8 along with a scatter-plot of the correlation for the 
individual pairs. Figure 9 displays plots for the variance and 
the spatial correlation for the KWAJ gauges for the JJA 
2001 and DJF 2002 seasons as illustration of the kind of results 
obtained. Since it includes only nine gauges, we display the 
correlation for the individual pairs without attempting to group 
them into spatial separation bins. Rough estimates of the point 
variance o\ and the corresponding values of the spectral cut-off 
A inferred from them using equation (34) are also included in 
the last two columns of Table 1. 

4. Discussion 

4.1. Radar Statistics 

[ 46 ] As already noted above, the complex behavior of the 
space-time statistics especially at the larger spatial scales 
limits the extent to which one can expect quantitative agree- 
ment between the radar data and our model. One limiting fac- 
tor is that the spectral model assumes space-time stationarity 
of the statistics at the outset. These assumptions greatly facil- 
itate the formulation of a statistical model in a parsimonious 
manner with a relatively small number of adjustable parame- 
ters. However, they also preclude application of the model to 
physical situations where local meteorological effects can, in 
general, lead to nonstationary statistics. The assumption of 
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spatial isotropy can be relaxed by introducing a form of r k in 
equation (6) that depends on both the magnitude and direc- 
tion of k. Such a model will necessarily contain more param- 
eters and will be more difficult to fit to data. 

[ 47 ] For both the MELB and the KWAJ data sets, the 
model is able to fit the spatial statistics fairly well at small 
spatial scales. The observed scatter in the plots of the spatial 
correlation function <b AA (s) at large s (Figures 3 and 4) in part 
reflects the fluctuation resulting from increasingly smaller 
sample sizes. However, the model has difficulty fitting the 
observed spatial correlation at larger separations where a 
small nonzero tail appears to linger. This is also reflected 
by the fact that the model prediction for the variance o\ as 
a function of box size L in general agrees well with the 
observed value for spatial scales up 32 km but falls short at 
the larger scales (see Figure 2). This is appropriate for an 
application of the model to problems that involve statistical 
comparisons between precipitation data from ground radar 
and a dense network of rain gauges located within the radar 
FOV. This is in contrast to BK96 where the fitting process 
was weighted so as to fit the larger scales better, since there, 
the intended application was to faithfully capture the statis- 
tics of large-area averages measured from satellites. 

[ 48 ] However, fitting the temporal statistics proved to be a 
rather challenging task. Quality of the fit for the temporal au- 
tocorrelation varies from one data set to another. Overall, the 
fit was better for the KWAJ data set than the MELB data set. 
The observed <&aa( r) often deviates from the model behavior 
at large spatial and temporal scales as evidenced by the pres- 
ence of heavy tails in Figures 5 and 6. Consequently, the fit 
had to be restricted to a limited range of space-time scales. 

[ 49 ] For the MELB radar, the complication arises in part 
due to the fact that it is located close to the Atlantic coast, 
and the 128 km box centered at the radar lies partly over land 
and partly over the open ocean. Thus, strictly speaking, the 
key assumption of spatiotemporally stationary statistics, 
which is at the heart of the theoretical foundation of the 
model, is violated here. This probably explains why the 
MELB data fits are in general the poorer of the two radar data 
sets. We have already noted that the lagged autocorrelation 
®aa(j) exhibits oscillatory tail behavior at large lags that be- 
come more pronounced for the larger spatial averaging scales. 
Our spectral model also predicts a similar behavior in the case 
/? > 1 (see Figure Bl) that is, however, too small to match the 
observations. This makes it rather difficult to find a set of tem- 
poral parameters that fit O aa (z) uniformly well over the entire 
range of scales. Although we are unable to fully explain this 
discrepancy, we suspect that it is partly attributable to the inter- 
action between a strong local diurnal cycle and the land-sea 
contrast in the coastal region. Because of these oscillations, it 
is difficult to obtain a stable estimate of the integral correlation 
time z A from the data that could be reliably compared with the 
model prediction, as was done in KB03 for radar rainfall data 
over the open oceans to fit the model in the f}=\ case. For the 
KWAJ radar, the purely oceanic environment makes the as- 
sumption of spatial homogeneity a reasonable one. However, 
deviation from temporally stationary behavior is fairly com- 
mon even in an oceanic environment. This generally manifests 
itself in the form of a time-dependent mean and variance. In 
this case different subsets of a season, say different months, 
will have unequal rain rate variance. Fitting the model to data 
that is inherently heteroscedastic forces one to restrict the 


space-time size of the sample. Thus, there is an inevitable com- 
promise: too large a sample may strongly violate the model as- 
sumption of space-time stationarity, but too small a sample 
may yield noisy and unreliable statistics. Partitioning the data 
into seasons is merely a practical way to cope with the situa- 
tion, the working assumption being that the statistics remain 
stationary in the course of a particular season. Seasonal depen- 
dence of the model parameters accounts for variation of the 
statistics from one season to another. This may not always 
work as the KWAJ DJF 2001 data illustrates, where isolated 
events cause substantial deviation from stationarity. Because 
of temporally nonstationary behavior, the statistics can defy 
description in terms of a spectral model like one described in 
this paper. The DJF 2002 results presented here were relatively 
well-behaved. Nonstationary statistics caused by local meteo- 
rological effects, such as land-sea contrast and orography, are 
rather difficult to quantify. A model like ours has no effective 
way of dealing with it and thus the prospect of a good agree- 
ment is somewhat tempered in many cases. 

[ 50 ] Clearly, the radar data does not currently have adequate 
spatial resolution to be able to reveal substantial departure 
from the power law singular nature of the spatial statistics 
predicted by the model. In fact, the singularity of o\ is weak 
enough to be still consistent with its observed growth rate at 
the smallest length scales probed by the radar. Consequently, 
there is no need for postulating the spectral cut-off A discussed 
in section 2.5 purely on the basis of the radar statistics alone. 
However, the need for a small distance cut-off becomes appar- 
ent when one considers the statistics of gauge data, which can 
access a broader range of space-time scales. 

4.2. Gauge Statistics 

[ 51 ] Wang et al [2008] recommend that while comparing 
the model prediction with the observed gauge statistics, one 
should exclude the values of time averaging scale T less than 
4-7 min so that the TB gauge estimates can be considered to 
be reasonably accurate. In Figure 7 we note that the model 
substantially underestimates (by a factor of about 1.5 to 2) 
the observed values of o\ obtained by averaging over all 
the available gauges in the MELB radar FOV. In Figure 8 
the spatial correlation W TI {p) for individual gauge pairs 
exhibits a rather large scatter for all the four 2001 seasons 
examined, perhaps indicating substantial spatial anisotropy 
of the rain statistics in the MELB radar area. The model repro- 
duces the binned average of the observed values of W TT (p) 
reasonably well for small separations. The anomalous 
behavior of the observed values appearing at longer distance 
scales is presumably caused by large-scale phenomena similar 
to those in the radar data, and like the latter, it is difficult to 
associate it with local meteorology. The large scatter in the 
values of the spatial correlation for the individual gauge pairs 
in the vicinity of the MELB radar (Figure 8) is similar to what 
was seen in [ Ciach et al ., 1997] for daily averages near 
Darwin, Australia. Agreement between the observed and the 
model-generated statistics for the gauges is considered 
adequate given the limitation that the gauges only sample the 
land portion of the 128 km box. For the KWAJ data 
(Figure 9), far fewer data points are available, but the situation 
appears to be generally similar. The observed values of o\ 
agree fairly well with the model prediction, but the extent of 
agreement is less clear for the spatial correlation W T1 {p). 
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[ 52 ] For averaging time scales Tless than about 30 min, the 
observed g\ does not increase as steeply as predicted by the 
spectral model, for which the asymptotic behavior (B26) 
suggests a power law divergence. Rather, it appears to flatten 
out to a finite value o\ as we consider smaller and smaller 
values of T. While this may in part be due to a possible distor- 
tion of the TB gauge statistics at small T, we suspect that this 
effect alone is not enough to account for the observed modi- 
fied behavior. In fact a similar flattening trend is noticeable 
also for the variance of the rain rate as measured by an optical 
disdrometer (courtesy Dr. Ali Tokay, not shown), which does 
not suffer from the kind of error that plagues the TB 
estimates. It appears reasonable to conclude that g\ actually 
tends to a finite value g\ as T 7 — ► 0 indicating the need of a 
short distance cut-off A. As expected, the values of A inferred 
from the rough estimates of g l listed in Table 1 are much 
smaller than the smallest radar spatial scale, 2 km, thus 
providing an a posteriori justification for neglecting it during 
parameter estimation from radar data. 


historical developments. For a noninteger order f (which 
can in general even be complex), the usual integral and deriv- 
ative operators are generalized into one-parameter families of 
operators and labeled by a real variable a , sometimes 
collectively referred to as fractional differintegral operators 
[Oldham and Spanier , 2006]. An integral of fractional-order 
a is defined by the formula 

J a t f(t) = - f\ a Auf(u)(t - (a > 0). (Al) 

In the case a = n > 0, a positive integer, it correctly reduces to 

Ct f Ml j* Un - 1 

= J a <*M 1 J a d «2 •••.!<, d U n f(u n ) 

= diif{u)(t - u) n ~\ 

which is the unique solution to the differential equation 


5. Summary and Conclusion 

[ 53 ] In this paper we have presented a spectral model of 
precipitation that describes the second moment statistics of 
rain rate from both radar and rain gauge data within a unified 
framework. The model is based on a linear stochastic differ- 
ential equation involving a fractional-order time derivative 
that generalizes a simpler predecessor based on an equation 
that was first order in time. The parameters of the model spec- 
trum are estimated by fitting the radar data. Our parameter es- 
timation can be thought of as being a “skewed fit” since the 
optimization procedure we adopted is designed to better esti- 
mate the smaller scale property of the precipitation field. 
With two adjustable spectral indices, the model fits the ob- 
served statistics of radar rainfall fairly well at smaller spatial 
and temporal scales. The model is tested by comparing its 
predictions against the observed statistics of rain rate data 
taken from a set of gauges located within the radar FOV. 
With the radar parameters, the model appears to reproduce 
the scale dependence of the gauge statistics reasonably well. 
The difficulties experienced by the model at larger scales vary 
from site to site and from one season to another. They are pre- 
sumably attributed to the statistically nonstationary nature of 
the data that is intrinsic to the rainfall phenomenon itself. 

[ 54 ] The spectral model provides a useful tool for estimation 
of the sampling error for various radar-gauge intercomparison 
scenarios [Bell and Kundu , 2003]. Radar at a GV site that is 
equipped with a dense rain gauge network can be conveniently 
calibrated using the model as a framework for data fitting. We 
plan to explore this problem further in a future paper. 


Appendix A: Fractional Calculus — A Brief Outline 

[ 55 ] The notion of a derivative or integral of fractional order 
has a long history that is almost as old as calculus itself going 
back to Leibniz, who mentioned a quantity that in modem ter- 
minology would be equivalent to derivative of order 1/2. The 
concept was first systematically formulated by Liouville with 
important subsequent contributions by Riemann, Gninwald, 
Krug, Sonine, Weyl, Riesz, and many other mathematicians. 
See, e.g., Samko et al. [1993] for a detailed account of the 


y (n \t ) =f{f),y{a) =y(a) =y"(a) = ... =/ n ~ l \a) = 0. 

The derivative operator J^ t is then defined through 
1 


J duf(u) R p o 

■tf/0= / v } ' ^' M) 

(<y) aI " “ 1 < Rey5 < "> n > °- 


(A2) 


[56] There is a reciprocal relationship between the frac- 

tional integral and derivative operators that is suggested by 
analogy with traditional calculus. It is usual practice to con- 
solidate the notation by writing J^ t = thus formally 

extending the fractional derivative operators to all orders in 
a way that clarifies the inverse relation between the two sets 
of operators introduced separately. 

[ 57 ] We now list some of the basic properties: 

1. a I“-/f(t) = = /rami 

2 . a m) =/(o, aim) = & m = f (n \t ), 

3. jrj{t) = 

The last relation indicates that the operators d/d t and a 
l)~ a do not commute. 

[ 58 ] Two cases can be distinguished. In the case when a is 
finite, one can choose a = 0 without loss of generality. The 
derivative operator 0 Z)f is the so-called Riemann-Liouville 
derivative , which acts on classes of functions defined on the 
half line [0, oo) and naturally appears in the context of an 
initial value problem like Brownian motion. The other 
case arises when one lets a — > — oo. The derivative _ 00 £>f , 
which we shall call the Liouville-Weyl derivative , emerges 
as the relevant fractional derivative operator when one is in- 
terested in classes of functions that reside on the entire real 
line (— 00 , oo). (Weyl’s name is attached to it because he 
introduced a class of fractional operators that coincides with 
_oo Z>f for the class of periodic functions). Now the initial 
condition refers to the infinitely distant past and the temporal 
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Figure Bl. Plot of the function h(rj) versus r\ for several 
values of the parameter /?. 


evolution is time translation invariant. It is the latter case 
that is of greater interest to us, since it allows one to describe 
temporally stationary stochastic processes. 

[ 59 ] The Liouville-Weyl derivative has the important 
property that under the operation of Fourier transform 
defined by the representation 

a(k, t ) = (27r) _1 ^ 2 j*_ 00 d co e~ lcot a{ k, co ), 


Plots of h{rj) are shown in Figure Bl for several values of 0 
between V 2 and 2. The case 0 = 1 is elementary: the function 
h(rj) is simply exp(— rf). When 0 < 1, h(rj) decreases mono- 
tonically and rapidly approaches zero for large r\. When 
0> 1, h{rj) exhibits damped oscillatory behavior around 
zero as rj increases. 

[ 62 ] Next we turn to the statistics of area-averaged rain 
rate. We start with the definition of the lagged rain rate 
covariance of two areas 


T^(s,r) = (l/T 4 )J d 2 x J d 2 x' c(s + x' - x,r), (B 3 ) 

A A' 

which can also be expressed in terms of a Fourier integral 
representation 


r^(s,r)=( 2 / 7 r)Jo d^J 0 dk 2 e zk * s sine 2 (k\L/2) sine 2 {k 2 L/2)c(k, r), 

(B 4 ) 

where sinc(x) = sin(x)/x and k = (k\ + kty ^ 2 . In the case of 
zero lag, i.e., r = 0, following the steps detailed in BK96, 
we obtain the spatial covariance function 

r^(>,0)=y 0 J_ 1 d4;J_ 1 df 2 (l-|f 1 |)(l - \i 2 \)C v {% + s/L\{L/L,)). 

(B 5 ) 

For zero separation (s = 0), we get the temporal 
autocovariance of an area A: 


the usual correspondence d/d t <=> — ico for ordinary de- 
rivatives holds for the fractional derivative as well; that is 

to say, if f{t) and F(co) are Fourier transforms of each other, 
then and (— icofFigo ) are Fourier transforms of 

each other as well. 


Appendix B: Derivation of Some Spectral 
Model Formulas 


[ 60 ] In this Appendix we describe some of the details of the 
analytical steps in the derivation of the various covariance 
statistics of precipitation data that are derived from our frac- 
tional dynamics model. Most of the steps are very similar to 
those in the 0 = 1 case of the present model described in 
BK96 and KB03 and we will be brief about them. 

[61] First, we examine the lagged covariance of the Fourier 
components of the precipitation field c(k , r) given by equa- 
tion (13), which is the basic building block of the various 
temporal statistics of interest. The function h{rj) is defined 
by the integral representation 


h(rj) 


' 2\ 1/2 M cos (fr) 

FJ g(P) J ° <f 2/? + 2 cosGfor/2)^ + 1 


(Bl) 


with h{ 0) = 1 . An explicit formula for the normalization factor 
gift) is obtained by setting rj = 0 and evaluating the f -integral. 
When — 2 < < 2, it has the explicit form 


Taa( 0,r) = (2 / 7r) J 0 d£ij 0 d k 2 sine 2 (k \L/ 2) sine 2 (k 2 L/ 2) c(k, z). 

(B6) 

Figure B2 illustrates the improvement in the fit to the lagged 
autocorrelation function 0^(0, r) over the/? = 1 case in BK96 
achieved by letting / 3 be an adjustable parameter. 

[ 63 ] The variance of area-averaged rain rate T^(0, 0) 
= o\ can be expressed in two alternative ways, as a spatial in- 
tegral obtained by setting the separation s = 0, in equation 
(B3) and as a Fourier integral by setting r = 0 in equation 
(B4). The first approach yields the formula 


^<^}= 4 y 0 G(v;L/Z 0 ), (B 7 ) 


where we have defined the function 


G(v;z) = J 0 dfi 0 df 2 (l - f,)(l - {2)Cy(z(^$+&y). 

(B 8 ) 

As pointed out in KB03, upon equating the two formulas in 
the case 0= 1, one obtains the useful identity: 

f°° f 00 2 71 

Jo dtfJo &k 2 I (k\ ,k 2 ;L/Lo, X) = r ^ ^ G( 2 ;T/T 0 ), (B 9 ) 

where, for brevity, we have defined 


yfhz eot(fln/2) 
0 sin(7r/y^) 5 


(B 2 ) 


which can be shown to be positive despite its appearance. 


I(k u k 2 ,L/L 0 i X) 


sine 2 (k\L/ 2L 0 ) sine 2 ( k 2 L / 2 L 0 ) 


[l + k\ + kI\ 


1+2 


(BIO) 
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Figure B2. Sample plots of the observed temporal autocor- 
relation ® aa (t) for a 128 km box (open circles) compared 
with its model prediction corresponding to the best fit values 
of P and zq (J3 = 1 .33, r 0 = 750min) for the same spatial scale 
(solid curve) and in the case ft=l with best fit value r 0 = 2460 
min (dashed curve). 


Note that v' = v in the “old” spectral model with /? = 1. From 
equations (16) and (B14) we have the general relation 


1 _ 1 +v 
~P ~ 1 T V ‘ 


(B15) 


From equation (B15) it follows immediately that v' > v and 
v'<vv correspond respectively to the cases P > 1 and P < 1 . 
Combining equation (B13) and the expression (B7) for a A , 
we get an explicit formula for z A : 


ft To r(l+v) G(1 + 2v';L/Lq) 
ZA ~\2g{fi)T(2 + 2v) G(v;L/Lo) 


(B16) 


[65] Next we outline the computation of the temporal sta- 
tistics g 2 t and 'i'riip) needed to describe the gauge data. 
The space-time covariance of rain rate averaged over two 
time intervals of length T lagged by time r measured at two 
points x and x f separated by a distance p, namely, 

c t 

r Tf {p,T)m((l/T)] 0 dtR' 

Ct Ct 

= (l /T 1 )] 0 d/J 0 dt'c{p,t - f + t) 


Ct 

(x, 0-(l/^)J 0 d^'^^xC + r )> 


[64] Next we evaluate the characteristic time scale z A 
defined by equation (23) as follows. In KB03 the following 
identity was proved: 


can be evaluated as follows. First, as usual, we convert the 
double integral into a single integral yielding the expression 


t a o\ m s/njlY^w = 0), (Bll) 

where t AA (co) is the Fourier transform of 1^(0, r). In 
view of the Fourier representation (B6), T AA (co) can be 
expressed in a similar form 

|*co Coo 

T 'aa(w) — (2 / tt) J 0 d^ij 0 d ^2 sine 2 (k\L/ 2) sine 2 (k 2 L/ 2) S(k, co), 

(B12) 

where S(k,co) is the model spectrum. Evaluating the right hand 
side of (B1 1) with the explicit form of the spectrum and using 
the relation (18) between the constants F 0 and y 0 , we get 

Taa{oj = 0) 27 ° r (| + v ) J q dKiJ dK 2 I(K U K 2 ‘,L/L 0 ,aP - 1). 

Kg(P) 

With the help of the identity (B9) this can be recast in the 
form 


Yaa(co = 0) 


4y 0 r 0 r(l + v) 
T(2 + 2v')g(P) 


G(1 + 2v';L/Lq). 


Then equation (Bll) yields 


T 

r rr'(P>6 = ( 1 / 7 ’)l dr '( 1 “ y) c (P’ t ' ~ t ) 

T 

= (l/r)|dr'(l: - T -^j[c(p,i' -t) +c(p,t' + t)\. 

Next we introduce the space-time Fourier transform repre- 
sentation of c(p,z) via equation (11). After performing the 
r'-integration explicitly and some further simplification, we 
obtain the formula 

|*GO C CO 

r TT ' (/>, z)=\/2/7rJ 0 d co cos (cot) sine 2 (coT/ 2) J 0 dk kJ$(kp)S(k, co). 

(B17) 

At zero lag this reduces to the spatial covariance of time-aver- 
aged rain rate: 

C co C 00 

Ttt(p,0) = a/2 / 7i J 0 d co sinc 2 (ft>r/2)J 0 dk kJo(kp)S(k,co). 

(B18) 

The variance o\ is obtained by specializing equation (B 1 8) to 
the casep = 0: 


T AG A 


2\/27ry 0 T 0 r(l + v) 
T(2 + 2v')g(P) 


G(1 + 2v’;L/Lq), 


(B13) 


c t 2 r = a/2/7tJ 0 d co sine 2 (coT / 2)) q dk kS(k, co). 


(B19) 


where we have now introduced a new parameter V through 
the relation 

ap = 2(1 +v'). (B14) 


[66] Finally, we investigate the asymptotic behavior of the 
area and time-averaged statistics by isolating the leading sin- 
gularity in the area/time integrals in which c(p , r) is replaced 
by its scaling approximation c (co) (p, r). It is convenient to ex- 
press the singular term in terms of the dimensionless variables 
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P*,t*,L*=L/L 0 and T* = T/z 0 . In the scaling limit r*,Z*— >0, 
the lagged autocovariance of area-averaged rain rate has the 
singularity structure 

1^(0, r) ~ y 0 L; 2 M<p i (r*/Z>,/?) , (B20) 

where the scaling function (p x (u\a , /?) is defined as 

Pi(«;«,A)=4| 0 lo dz i dz 2(l-zi)(l-z2)(^+z^) |v| (9^(z^+4) a/2 ;a,{i). 

(B21) 

In particular, as was shown in KB03, the variance = 
T^(0,0) has the limiting behavior 

g 2 a ~A + BL- 2 |v| (B22) 

plus terms that vanish as L* — > 0. This result continues to hold 
in the more general /? f 1 case, where A and 5 are constants 
given by 

A = ^o r (-M)> 

5 = 2 1+2 M 70 r(|v|)J 0 dfJ 0 df 2 (l -£,)(1 -f 2 )(f 2 +f 2 ) v . 

(B23) 

[67] The integral correlation timescale given by equation 
(B12) has the power law dependence ta~L 2 \ v \ as Z*— >0. 
Similarly, in the scaling limit /?*, Z* — ► 0 the spatial covariance 
of time-averaged rain rate has the singularity structure 

TttCo, o) ~ 7o/U l 'V 2 (r./pj; «,/?) , (B24) 

where we have introduced a second scaling function 

(p 2 (s;a,P) = 2 # f 0 d^(l - z)<p(zs;a,fl). (B25) 

The asymptotic behavior of the variance of time-averaged 
rain rate o\ = T rr (0, 0) for small Z* is inferred from equations 
(B24) and (B25) by making use of the limit (32): 

4~(1 - M/«) _1 (l -2\v\/a)- l y 0 KT; 2 ^ a . (B26) 

[68] Acknowledgments. We thank David Marks of Science Systems 
Applications International, Inc. and NASA/GSFC for in-depth discussions on 
various details regarding the TRMM GV data used in this study and for provid- 
ing us with Figure 1 . We are also grateful to Bikas K. Sinha of Department of 
Mathematics and Statistics, UMBC for his careful reading of the paper and 
Anindya Roy also of Department of Mathematics and Statistics, UMBC for 
many help fill discussions. Comments and constructive criticisms from the 
associate editor and two other anonymous reviewers on various aspects of 
the work, especially the parameter estimation method, greatly helped to 
improve the presentation of the paper. The research conducted in this paper 
was supported by the GPM (Global Precipitation Measurement) GV group 
as a part of the NASA Precipitation Measurement Missions (PMM) program. 
One of us (J.E.T.) gratefully acknowledges financial support from a 
JCET Fellowship. 


References 

Bell, T. L. (1987), A space-time stochastic model of rainfall for satellite 
remote sensing studies, J. Geophys. Res., 92, 9631-9643, doi:10.1029/ 
JD092iD08p0963 1 . 

Bell, T. L., and P. K. Kundu (1996), A study of sampling error in satellite 
rainfall estimates using optimal averaging of data and a stochastic model, 
J. Climate, 9, 1251-1268, doi: 10.11 75/1 520-0442(1 996)009< 1251: 
ASOTE>2.0CO;2. 

Bell, T. L., and P. K. Kundu (2003), Comparing satellite rainfall estimates 
with rain gauge data: Optimal strategies suggested by a spectral model, 
J. Geophys. Res., 108{D3), 4121, doi:10.1029//2002JD002641. 

Beran, J. (1994), Statistics of Long-Memory Processes, Monographs on 
Statistics and Probability 61, Chapman and Hall, New York. 

Ciach, G. J., W. F. Krajewski, E. N. Anagnostou, M. L. Baeck, J. A. Smith, 
J. R. McCollum, and A. Kruger (1997), Radar rainfall estimation for 
ground validation studies of the Tropical Rainfall Measuring Mission, 
J. Appl. Meteor ol. , 36, 135-141, doi: 10.1 175/1520-0450-36.6.735. 

Efron, B. (1981), Nonparametric standard errors and confidence intervals 
(with discussion), Can. J. Stat., 9, 139-172. 

Frisch, U. (1995), Turbulence : The Legacy ofA.N. Kolmogorov, Cambridge 
Univ. Press, Cambridge, U.K. 

Habib, E., W. F. Krajewski, and A. Kruger (2001), Sampling errors of 
tipping-bucket rain gauge measurements, J. Hydrol. Engng., 6(2), 159-166, 
doi:10.1061/(ASCE)1084-0699(2001)6:2(159). 

Hawkins, D. L. (1989), Using U-statistics to derive the asymptotic distribu- 
tion of Fisher’s Z statistic, The American Statistician, 43(4), 235-237, 
doi: 10.2307/2685369. 

Jenkins, G. M., and D. G. Watts (1968), Spectral Analysis and its Applications, 
Holden-Day, San Francisco, CA. 

Kedem, B., and L. S. Chiu (1987), Are rain rate processes self-similar?, Water 
Resour. Res., 23(10), 1816-1618, doi:10.1073/pnas.84.4.901. 

Kundu, P. K., and T. L. Bell (2003), A stochastic model of space-time vari- 
ability of mesoscale rainfall: Statistics of spatial averages, Water Resour. 
Res., 39(12), 1328 doi: 10.1 029/2002 WR00 1802. 

Kundu, P. K., and T. L. Bell (2006), Space-time scaling behavior of rain sta- 
tistics in a stochastic fractional diffusion model, J. Hydrol., 322, 49-58, 
doi: 1 0. 1 0 1 6/j .jhydrol.2005 .02 .03 1 . 

Kundu, P. K., and R. K. Siddani (2007), A new class of probability distributions 
for describing the spatial statistics of area-averaged rainfall, J. Geophys. Res., 
112, D18113, doi:10.1029/2006JD008042. 

Lagarias, J. C., J. A. Reeds, M. H. Wright, and P. E. Wright (1998), 
Convergence properties of the Nelder-Mead Simplex Method in low dimen- 
sions, SLAM J. Optim., 9, 112-147, doi:10.1 137/S1052623496303470. 

Matem, B. (1960), Spatial variation: Stochastic models and their application 
to some problems in forest surveys and other sampling investigations, 
Medd.fran Statens Skogsforskningsinst., 49, 5, 1-144. 

Matem, B. (1986), Spatial Variation: Lecture Notes in Statistics, vol. 36, 
Second ed., 151 pp. Springer-Verlag, Berlin and New York. 

Miller, K. S., andB. Ross (1993 ), An Introduction to the Fractional Calculus 
and Fractional Differential Equations, John Wiley, New York. 

Nelder, J. A., and R. Mead (1965), A simplex method for function minimi- 
zation, The Computer Journal, 7, 308-313. 

North, G. R., and S. Nakamoto (1989), Formalism for comparing rain esti- 
mation designs, J. Atmos. Oceanic Tech., 6, 985-992, doi: 10.1 175/1520- 
0426(1989) 006<0985:FFCRED>2.0.CO;2. 

Oldham, K. B., and J. Spanier (2006), The Fractional Calculus, Dover 
Publications, Mineola, New York. 

Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery (1992), 
Numerical Recipes in FORTRAN: The Art of Scientific Computing, 
2 nd edition, Cambridge University Press, Cambridge, New York. 

Samko, S. G., A. A. Kilbas, and O. I. Marichev (1993), Fractional Integrals 
and Derivatives, Gordon and Breach Science Publishers, New York. 

Stom, R., and K. Price (1997), Differential evolution — A simple and effi- 
cient heuristic for global optimization over continuous spaces, J. Global 
Optimization, 11, 341-359, doi:10.1023/A:1008202821328. 

Wang, J., B. L. Fisher, and D. B. Wolff (2008), Estimating rain rates from 
tipping-bucket rain gauge measurements, J. Atmos. Oceanic Tech., 25, 
43, doi:10.1 175/2007JTECHA895.1. 

West, B. J., M. Bologna, and P. Grigolini (2003), Physics of Fractal 
Operators, Institute of Nonlinear Science, Springer, New York. 

Wolff, D. B., D. A. Marks, E. Amitai, D. S. Silberstein, B. L. Fisher, 
A. Tokay, J. Wang, and J. L. Pippitt (2005), Ground validation for the 
Tropical Rainfall Measuring Mission (TRMM), J. Atmos. Oceanic Tech., 
22, 365-379, doi: 10.1175/JTECH 1700.1. 


10,295 


